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Abstract 



Based on a re- formulation of the classical explanation of quantum mechanical Gaussian dispersion 
(Grossing et al. 2010 [.]) as well as interference of two Gaussians (Grossing et at 2012 [ ]), we 



present a new and more practical way of their simulation. The quantum mechanical "decay of 



the wave packet" can be described by anomalous sub-quantum diffusion with a specific diffusivity 

in 

qq varying in time due to a particle's changing thermal environment. In a simulation of the double-slit 

OO 

experiment with different slit widths, the phase with this new approach can be implemented as a 

local quantity. We describe the conditions of the diffusivity and, by connecting to wave mechanics, 

i—i we compute the exact quantum mechanical intensity distributions, as well as the corresponding 

trajectory distributions according to the velocity field of two Gaussian wave packets emerging from 

a double-slit. We also calculate probability density current distributions, including situations where 

a 

phase shifters affect a single slit's current, and provide computer simulations thereof. 
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1. INTRODUCTION 

In reference [ ] we presented a classical model for the explanation of quantum mechanical 
dispersion of a free Gaussian wave packet. In accordance with the classical model, we shall 
now relate it more directly to a "double solution" analogy gleaned from [ ]. For, as is 
shown, e.g., in [4, 5], one can construct various forms of classical analogies to quantum 
mechanical Gaussian dispersion. Originally, the expression of a "double solution" refers to 
an early idea of [6] to model quantum behavior by a two- fold process, i.e., by a the movement 
of a hypothetical point-like "singularity solution" of the Schrodinger equation, and by the 
evolution of the usual wave function that would provide the empirically confirmed statistical 
predictions. Recently, [3] used this ansatz to describe the behaviors of their "bouncer"- (or 
"walker"-) droplets: On an individual level, one observes particles surrounded by circular 
waves they emit through the phase-coupling with an oscillating bath, which provides, on 
a statistical level, the emergent outcome in close analogy to quantum mechanical behavior 
(like, e.g., diffraction or double-slit interference). The simulation of interference in the 
double-slit experiment was in [ ] easily achieved by assuming the simple case where the two 
slits have equal aperture. Instead, in this paper, we discard that simplification and show 
in a more detailed analysis that one can i) find a formulation applicable to independently 
variable slit widths, and ii) provide computer simulations thereof. 

2. FROM CLASSICAL PHASE-SPACE DISTRIBUTIONS TO QUANTUM ME- 
CHANICAL DISPERSION 

In the context of the double solution idea, which is related to correlations on a statistical 
level between individual uncorrelated particle positions x and momenta p, respectively, we 
consider the free Liouville equation 

dt ^-^ m dxi ^ dxi dp* 

4 = 1 4 = 1 

with potential V and mass m. For simplicity, we restrict ourselves to the 1-dimensional 
space coordinate x further on. Liouville's equation (2.1) implies the spatial conservation law 
and has the property that precise knowledge of initial conditions is not lost in the course of 
time. That is, it provides a phase-space distribution / (x,p,t) that shows the emergence of 



correlations between x and p from an initially uncorrected product function of non-spreading 
( "classical" ) Gaussian position distributions as well as momentum distributions, 

/o( "' p) = 2^ eXP {-^f} eXP {-^f}' (2 ' 2) 

where <To is the initial space deviation, i.e., ctq = a(t = 0), and tc q := mu is the momentum 
deviation. Then the phase-space distribution reads as 



1 I (x — pt/m) | f p 

exp < k > exp < ., 

2naomuo 2ctq { 2m 2 UQ 



/(»,M) = ^»P -^3^ «P -i ■ < 23 > 



The above-mentioned correlations between x and p emerge when one considers the proba- 
bility density in x-space, which is given by the integral 

p{x ' t) = I ,dp= ^k exp {-^}- {2A) 

with the variance at time t given by 

a 2 = al + u 2 t 2 . (2.5) 

In other words, the distribution (2.4) describing a spreading Gaussian is obtained from a 
continuous set of classical, non-spreading, Gaussian position distributions whose momenta 
also have a non-spreading Gaussian distributions. One thus obtains the exact quantum 
mechanical dispersion formula for a Gaussian, as we have obtained also previously from our 
classical ansatz by relating different kinetic energy terms in our diffusion model [ ] . For 
confirmation with respect to that model we use the Einstein relation 



with the reduced Planck constant h = h/(2ir) and m being the particle's mass, and we note 

JL 

dx J 



that with (2.4), VP = -k-P = —\P and the usual definition of the "osmotic" velocity u 



one obtains 

VP xD 

u = u{ X) t) = -D^- = X —. (2.7) 

P <7 Z 

For the average initial value we find 



Uq := u 



= u(a o ,0) = — , (2.8) 



which turns out to be the same as the initial velocity at position x = <Jo at t = 0. This is 
a characteristic value for Gaussians, which we simply called "initial velocity" in our recent 



papers. It corresponds exactly to the velocity u at starting time t = at the trajectory that 
has distance £(0) = ±cr from the maximum of the Gaussian (Fig. 3.1). With Eq. (2.8) one 
can rewrite Eq. (2.5) in the more familiar form 

- 2 = -o 2 (l + if ) ■ P-9> 

Note also that by using the Einstein relation (2.6) the norm in (2.3) becomes the invariant 

expression 

112 . 
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reflecting the "exact uncertainty relation" [7]. 

3. SPREADING OF THE WAVE PACKET DUE TO A PATH EXCITATION 
FIELD 

We note that u/uq is a spreading ratio for the wave packet independent of x. This functional 
relationship is thus not only valid for the particular point £(£) = er(£), but for all x of the 
Gaussian. Therefore, one can generalize (2.9) for all x, i.e., 

a a I DH 2 

£(*)=£(0)- where ±-=Jl + —L. (3.1) 

o"o o"o y 0o 

In other words, one derives also the time-invariant ratio for the spreading 

m = m = const . (3 . 2) 

In Fig. 3.1 the spreading according to Eq. (3.1) is sketched. We can now try to implement our 
previous assumption that the "bouncer" particle is phase locked with its nonlocal diffusion 
wave field such that the Gaussian describing the diffusion process has long undulatory tails 
representing the locking in with the undulations of the zero-point field. In other words, 
we can now re-consider our classical simulations of Gaussian dispersion and double slit 
interference [2], respectively, by constructing from (2.4) a description of our "path excitation 
field" via the introduction of the amplitude R as product of a Gaussian (at rest in the re- 
direction) and a plane wave (in the y-direction), 

R (x, t) = (27nx 2 ) ~ 1/4 exp | - £-^ | cos (k y y) . (3.3) 




Figure 3.1: Bohm-type trajectories for a quantum particle with initial Gaussian 
distribution exhibiting the characteristics of ballistic diffusion 



The product (3.3) is factorisable for all t into x- and y-dependent functions, due to the 
motion in the ^/-direction by y (t) = hkt/m. According to our principle of path excitation 
[2], we deal with a single, classical particle of velocity v = p/m following the propagations of 
waves of equal amplitude R comprising a wave-like thermal bath that permanently provides 
some momentum fluctuations Sp. 

The latter are reflected in Eq. (2.9) via the r.m.s. deviation a(t) from the classical path. In 
other words, one has to do with a wave packet with an overall uniform velocity v , where the 
position x = vt moves like a free classical particle, as indicated in Fig. 3.1. As the packet 
spreads according to Eq. (2.9), £(£) describes the result of the motion along a trajectory 
of a point of this packet that was initially at £(0). The smaller the initial value of |£(0)|, 
i.e., the distance from x of the center point of the packet, the slower said spreading takes 
place. In our model picture, this is easy to understand: For a particle exactly at the center 
of the packet, x to t = ^o •"• £(0) = , the momentum contributions from the "heated up" 
environment on average cancel each other for symmetry reasons. However, the further off 
a particle is from that center, the stronger this symmetry will be broken, i.e., leading to a 
position-dependent net acceleration or deceleration, respectively, or, in effect, to the "decay 
of the wave packet" . The actual decay of the wave packet starts, roughly spoken, at a time 



tk, indicated by a "kink" in Fig. 3.1 which is due to the squared time-behavior in Eq. (2.9). 
From Fig. 3.1 we find x tot = x + v(t)t + f (t) and f (t) = f (0) + u(t)t. Without loss of 
generality we set v = const, and x = further on. With the use of Eq. (3.1) we obtain 



a \ v 2 t 2 

x tot {t) = vt+^t) = vt+m- = vt+m\h + ^r- (3.4) 

In our model picture, x tot is the position of the "smoothed out" trajectories, i.e., those 
averaged over a very large number of Brownian motions. 

Moreover, one can now also calculate the average total velocity field of a Gaussian wave 
packet as 

vUt) = %^ = v + £(0) , U ^l = , (3.5) 

which describes the velocity field v tot of a point along a trajectory (i.e, the residue of the 
"path excitation field" to be explicated further below). 

Finally, we derive the average total acceleration field of a Gaussian wave packet as 

q,o,W = %^ = «0) . <K , (3.6) 

describing the acceleration of a point along the trajectory at time t. Eqs. (3.4) to (3.6) allow 
us to calculate the quantities along a trajectory only from a given starting point, indicated 
by£(0). 

Actually, however, we are interested in the dynamics at any given position (x, t) directly. 
Using 

£(t) =x-vt (3.7) 

and Eq. (3.1) we rewrite 



4«»- * = (3-8) 

1 + U 2 Q t 2 /(T 2 



which leads to the generalized fields, 



^tot(^) t) - 


= X, 


Vtot(x,t) : 


--^^K u Sk^ + ^-<l 


a tot(^) t) '■ 


=«" (1+ ^r {x - vt) ^°' 



(3.9) 
(3.10) 
(3.11) 



which will be used in the simulations later on. 



4. THE DERIVATION OF D t 

We derive a solution for a diffusion equation with a time-dependent diffusion coefficient 
kt a for a generalized diffusion equation (cf. [ ]) 

dP d 2 P 

Here, t and fe denote the time and a constant factor, respectively. Inserting P(x, t) of 
Eq. (2.4) as a solution into Eq. (4.1) yields 

and, after integration, 

- = A;- + * . (4.3) 

2 a+1 2 v ; 

Substitution of (2.9) into (4.3) yields cq = a 2 and 

. 2t a+1 DH 2 , s 

(4.4) 



a + I a 2 
which can only be fulfilled by a = 1, so that 



D 2 

k = — . (4.5) 



o"6 



The time-dependent diffusion coefficient D t is with (2.8) identified as 



'0 

Finally, Eq. (4.1) reads as 

dP D 2 t d 2 P 



A: =5' = " Si = 5^ f - (46) 



dt o\ dx 2 



(4.7) 



and turns out to be a ballistic diffusion equation, defined by a = 1, as the special case of an 
anomalous diffusion where the diffusion coefficient D t grows linearly with time t. 

Essentially, the "decay of the wave packet" thus simply results from sub-quantum diffusion 
with a diffusivity varying in time due to the particle's changing thermal environment: as 



the heat initially concentrated in a narrow spatial domain gets gradually dispersed, so must 
the diffusivity of the medium change accordingly. 

Now we look at the time t^ of the kink (Fig. 3.1). The wave packet begins to spread 
differently at the kink, which is, according to Eq. (3.1), obviously at that time t = t\ where 
the influence of the right term is equal to the left term under the square root and hence 
DH 2 = g$ (i.e., o = y/2a ). Then we find with (4.6) that 

A = r D - ( 4 - 8 ) 

As one can see, t = t^ is the time when D t = D. Note that the diffusivity D is constant 
for all times t and has to be distinguished from the diffusion coefficient D t . In a different 
approach, one could also start out with the "exact" uncertainty relation, mu^t^ = H/2, with 
Mo = D/gq. This again leads to D t = DH/g 2 . = u^t = tjt^D. 

We recall Boltzmann's relation AQ = 2uqSS [9, 10] between the heat applied to an 
oscillating system and a change in the action function SS = ^5 J r E kin dt, respectively, 
providing 

VQ = 2oo V(5S) . (4.9) 

Here, oS relates to the momentum fluctuation via 

V{8S) = 8p=:mu=-~-, (4.10) 

and therefore, with P = P oe - s ^ kT ° and AQ = kT = hco, 

mu=^. (4.11) 

2u 

Using the initial velocity (2.8) together with Eq. (4.8) we find 

D t = u 2 t = \mu\W = 1^-tD = \ [SS(t) - 5S(0)} D = ^§^D. (4.12) 

lb 1 1 —, 1 1 I lb I bUJ 

Actually, 5S(0) = 0, since there are no initial fluctuations. Substitution of (4.8) into (4.12) 
leads then to 

Dt = l k D = ^D = -\n^D = -D[\nP(t) -lnP(O)]. (4.13) 

One can also derive a condition that does not require to know the diffusion coefficient at 
t = 0, 

AD = A(t 2 ) - A(*i) = ~D [lnP(t 2 ) - lnP(ti)] , (4.14) 



by choosing two arbitrary time steps t\ and £2 as suggested in Fig. 3.1. 
From condition (4.6) one can immediately see that 

dD t D 2 



dt a 2 



const. (4.15) 



Thus, one can also rewrite Eq. (4.14) as 



D 2 

AD = D t (t 1 ) + —(t 2 -t 1 ), (4.16) 



which is only valid for equal slit widths X01 = £02 ari< ^ thus °"i = °"2- In order to compute 
the distribution of P(x,t), one starts with Eq. (4.16) and takes the local properties of the 
diffusivity into account. For given times t± and £2 = ti + At one obtains with (2.8), 

A(Mi) = -^hiP(Mi), (4.17) 

D t {x,t 2 ) = -DhiP{x 1 t 2 ) + D t {x,t 1 ), (4.18) 

thereby constituting a rule to numerically compute the distribution P(x,t). 

5. FINITE DIFFERENCE SCHEME 

Starting with the ballistic diffusion equation (4.7) with time-dependent diffusivity D t we 
use an explicit finite difference forward scheme (cf. [11]), 

¥L^±(P[ x ,t + l]-P[x,i\), (5.1) 

^ -> -L ( P [ x + M ] - 2 P[x, t] + P[x - 1, t]) , (5.2) 

with 1-dimensional cells. In case D t is independent of x, the complete equation after re- 
ordering leads to 



D[t + l]At 
~Ax~ 2 



P[x,t + 1] =P[x,t}+ [ A 2 J — {P[x+l,t]-2P[x,t]+P[x-l,t}} (5.3) 



with space x and time t, and initial Gaussian distribution P(x, 0) with standard deviation 
a . 

As can be seen, calculation of a cell's value at time t only depends on cell values at the 
previous time. The time-dependence of the diffusion coefficient can also be calculated with- 
out any knowledge of neighboring cells. The diffusion coefficient represents the underlying 
physics of the current cell and is calculated for the evaluated time step t + 1. 



The stability condition for the scheme of Eq. (5.3) is that 

A A* 



< \ (5-4) 



L\X (Tn 

to < -^r ■ ( 5 - 5 ) 



Ax 2 

be satisfied for all values of the cells [x,t] in the domain of computation. The general 

procedure is that one considers each of the frozen coefficient problems arising from the 

scheme. The frozen coefficient problems are the constant coefficient problems obtained by 

fixing the coefficients at their values attained at each point in the domain of the computation 

(cf. [12]). Substituting Eq. (4.6) into (5.4) leads to 

Ax 2 aj 
~2Dfi 

This shows that the finite difference scheme (5.3) is suitable to solve the ballistic diffusion 

equation (4.7) as long as the spreading is not too big. Beyond that, the computations are 

no longer economically practical due to the necessarily enormous number of cells. Then, one 

has to replace the explicit scheme (5.3) by an implicit scheme, for example, which has less 

stringent restrictions on the stability conditions, but needs linear equation solvers instead. 

For our simulations we employed the explicit scheme introduced above as well as implicit 

schemes with an open source software for numerical computation, Scilab [ 3], on a standard 

personal computer. 

6. THE CONNECTION TO WAVE MECHANICS: THE DOUBLE-SLIT EXPER- 
IMENTS WITH DIFFERENT SLIT WIDTHS 

For a more generalized picture, we now take a closer look at the double-slit experiment. 
Consider a scenario as shown in Fig. 6.1 with two slits of different widths. We assume the 
initial Gaussians passing through a slit have a standard deviation value matching the slit 
width, e.g., o"oi = 0"o and cr 2 = 0"o/2, respectively, with o- Qi then being also the width of 
slit i. The resulting Bohm-type trajectories of the two decaying Gaussians are sketched in 
Fig. 6.1 with red lines. Thus 

t k2 =W4 (6.1) 

while the spreading is doubled (compare with the grayed out spreading of slit 2 for the case 
of o"o for both slits). According to Eq. (4.6), the diffusion coefficients of the two slits yield 

DH DH 

A,iW = — ^A, 2 (t) = -r, t>0. (6.2) 
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Figure 6.1: Sketch of a double-slit with two different widths and Bohm-type trajectories 
(and same-widths scenario indicated by gray lines) 



The advantage of Eq. (4.13) lies in it's local fit due to its dependence on P(x,t) instead 
on a(t), since the latter is just a global statistical value of too less local relevance. For 
the general case, we have to deal with a diffusion coefficient D t (x,t) further on. The time- 
dependent diffusion equation reads then 

OP d / . dP\ , , 

We have now all the tools necessary to consider the inclusion of wave mechanics in our 
model. We define the phase as 

tp = S/h (6.4) 

with the general action S. Identifying t> tot of Eq. (3.10) with 



Utot 



m 



we find for the action 



Edt 



m 



v + 



u 2 t 



Z(t)\ dx- / Edt 



(6.5) 



(6.6) 



al + ult* 

with E being the system's total energy. As v does not depend on x we can solve the first 
integral, and for the conservative case also the second integral, providing with (3.2) 



S = mvx + 



rnu 



a(t) 



t — Et = mvx + 



muz 






t-Et. 



(6.7) 
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Here, the action S along a trajectory is given by the sum of the usual momentum-related 
term and a term depending on the kinetic energy, or kinetic temperature, respectively, of the 
"heated up" environment, weighted by a factor that solely depends on a particular trajectory 
indicated by the initial location £(0) in the Gaussian. 

Finally, we rewrite the phase with the help of Eqs. (6.4) and (3.7) as 



<P 



1 
h 



mvx + 



muz 



O"0 



t- Et 



1 
h 



mvx + 



raur 



x — v t 



t-Et 



(6.8) 



The expression containing £(0) indicates a phase p along a trajectory, while the r.h.s. sticks 
to our coordinate system and is thus the better choice to do interference calculations. 

Instead of following just one Gaussian, we extend our simulation scheme to include two 
possible paths of a particle which eventually cross each other. For this, we use two Gaussians 
approaching each other. Following our earlier approach in [ ] we simulate a double-slit 
experiment by independent numerical computation of two Gaussian wave packets with total 
distribution given by 

P tot :=Pi + P 2 + 2 y^\P 2 cos (p 12 . (6.9) 

Since each Gaussian has its own phase (6.8) we are free to add a phase shifter A(p for one 
of the slits of the two-slit experiment, say slit 1, which modifies the phase to 

Si 



<Pi 



h 



+ Ay? 



(6.10) 



and yields for the phase difference 



<Pu = ¥2 - V?i 



m 
+ 



v 2 {x-x m ) -v x {x-x m ) 



mt 
2h 



ul 2 (x - X 02 - V 2 t) 2 V? m (x - Xqi - Vit) 



21 



- Ay?. 



(6.11) 



al(t) a\{t) 

The two slits at positions 201 an d £02 have different slit widths and hence different param- 
eters, coi, 01, tool an d 002, o~2, to 2, respectively, as illustrated by the red trajectories in 
Fig. 6.1 for the example of 001 = 00 and 002 = 0o/2, respectively. One can observe several 
characteristics of the averaged particle trajectories, which, just because of the averaging, are 
identical with the Bohmian trajectories. As one can see, the phase difference (6.11) is at 
any time defined for the whole domain, and hence <£> 12 is intrinsically nonlocal. 
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Finally, we recall our derivation of the total average density current [2, 14], i.e., the most 
general expression (including weights Pi) for our "path excitation field" , 

Jtot = P1V1 + P2V2 + \fpyPi {v\ + v 2 ) cos ifi2 + \fP~1P~2 (ui - u 2 ) sin ip 12 , (6.12) 

where 

Vtat = ~ET- ' ( 6 - 13 ) 

*tot 

with osmotic velocities Ui of Eq. (2.7) and total velocities Vi of Eq. (6.5) applied to both 
slits, 1 and 2, and with the phases (6.11). Note that the last term on the r.h.s in Eq. (6.12) 
is termed "entangling current" J e by us [ ], which is of a genuinely "quantum" nature in 
that the velocities Ui are generally entangled with the velocities t>;. 

7. SIMULATION RESULTS 

In Figs. 7.1 to 7.3, the graphical results of a classical computer simulation of the interfer- 
ence pattern in double-slit experiments are shown, including the trajectories. In Fig. 7.1a 
the maximum of the intensity is distributed along the symmetry line exactly in the middle 
between the two slits. [2] 

In the examplary figures, trajectories according to Eq. (6.9) for the two Gaussian slits are 
shown. The interference hyperbolas for the maxima characterize the regions where the phase 
difference if = 2mr, and those with the minima lie at ip = (2n + l)n, n = 0, ±1, ±2, . . . Note 
in particular the "kinks" of trajectories moving from the center-oriented side of one relative 
maximum to cross over to join more central (relative) maxima. In our classical explanation 
of double slit interference, a detailed "micro- causal" account of the corresponding kinematics 
can be given. The trajectories are in full accordance with those obtained from the Bohmian 
approach, as can be seen by comparison with [4, 15-17], for example. 

We use the same double-slit arrangements in Figs. 7.2 and 7.3, but include a phase shifter 
affecting the current from slit 1, as sketched by the dashed red lines or rectangles on the 
left hand side, respectively. Even though the total applied phase shift is either 37r or 5ir in 
Figs. 7.2a and 7.3a, respectively, one recognizes the effective phase difference of A(p mo d = Atp 
mod 27r = 7r in each case, which eventually results in equal shifts of the interference fringes. 
By comparing with Fig. 7.1a we now observe a minimum of the resulting distribution along 
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slit 1 slit 2 

(a) equal slit widths 



slit 1 



slit 2 



(b) cr i = 2<T ( 
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Figure 7.1: Classical computer simulation of the interference pattern with different slit 

widths: intensity distribution with increasing intensity from white through yellow and 

orange, with trajectories (red) for two Gaussian slits (v X) i = v Xt2 = 0) 



the central symmetry line, in full accordance with the Aharonov-Bohm effect, independently 
of the times t\ and t% during which the shift has been applied. 

To bring out the shifting of the interference fringes more clearly, we apply in Fig. 7.3 
the phase shift in between the indicated times t\ and t%, respectively. Note that the phase 
shift applied only to a single slit's current at a time t when the decaying Gaussians are 
already overlapping (Fig. 7.3) is shown here for didactic reasons only. In this highly idealized 
scenario, then, one can see an illustration of the immediate effectiveness of Atp over the whole 
domain according to Eq. (6.11), i.e., of the nonlocality of the relative phase. 

To conclude, we have in this paper provided a detailed description of the velocity fields 
involved in the analytical calculations as well as the computer simulations illustrating Gaus- 
sian dispersion and interference at a double slit. We have arrived at an expression for the 
local value of the phase, Eq. (6.8), which made it possible also to extend our previous model 
to slit systems with independently variable slit widths. With the computer simulations of 
the latter, the nonlocal nature of the relative phase can be clearly demonstrated. 
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slit 1 slit 2 

(a) Probability density P 
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(b) Phase shift 
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slit 1 slit 2 

(c) Entangling current J e 



Figure 7.2: Classical computer simulation as in Fig. 7.1a, but with additional phase shift 
Aip = 3n accumulated during the time interval between ti and t 2 at slit 1 




slit 1 slit 2 

(a) Probability density P 



5n A f 

(b) Phase shift 



slit 1 slit 2 

(c) Entangling current J e 



Figure 7.3: Classical computer simulation as in Fig. 7.2, but with different times tj and 

with accumulated additional phase A(f = 5tt. This results in the same distributions of P 

and J e for times t > t 2 and shows the effect of the shifting of the interference fringes more 

clearly than Fig. 7.2 
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